Magnetic Anomaly Maps
This file is best viewed in a Pluto notebook. To run it this way, from the MagNav.jl directory, do:
julia> using Pluto
julia> Pluto.run() # select & open notebookThis is a reactive notebook, so feel free to change any parameters.
Import packages and DataFrames
The DataFrames listed below provide useful information about the flight data (collected by Sander Geophysics Ltd. (SGL) in 2020) & magnetic anomaly maps.
| Dataframe | Description |
|---|---|
df_map | map files relevant for SGL flights |
df_cal | SGL calibration flight lines |
df_flight | SGL flight files |
df_all | all flight lines |
df_nav | all navigation-capable flight lines |
df_event | pilot-recorded in-flight events |
begin
cd(@__DIR__)
# uncomment line below to use local MagNav.jl (downloaded folder)
# using Pkg; Pkg.activate("../"); Pkg.instantiate()
using MagNav
using CSV, DataFrames
using Plots: plot, plot!
using Random: seed!
using Statistics: mean, median, std
seed!(33); # for reproducibility
include("dataframes_setup.jl"); # setup DataFrames
end;
Load Perth map
This is the Perth map (at 800 m) as provided by Sander Geophysics Ltd.
begin
map_gxf = MagNav.ottawa_area_maps_gxf()*"/Perth_Mag.gxf"
p_mapS_800 = map_gxf2h5(map_gxf,800)
p1 = plot_map(p_mapS_800)
end
Display Perth map in Google Earth
Display the Perth map in Google Earth by uncommenting below to generate a KMZ file, then open in Google Earth.
# map2kmz(p_mapS_800,"Perth")
Overlay Perth mini-survey
Overlaid on the Perth map are most of the flight lines used to generate the map, which were collected during Flight 1004 (see readme) & Flight 1005 (see readme). The full list of SGL flights is in df_flight.
begin
xyz_1004 = get_XYZ(:Flt1004,df_flight;silent=true) # load flight data
xyz_1005 = get_XYZ(:Flt1005,df_flight;silent=true) # load flight data
lines = [4019,4018,4017,4016,4015,4012,4011,4010,4009,4008,4007,4004,
4003,4002,4001,421,419,417,415,413,411,409,408,407,405,403,401]
for line in lines
xyz = line in xyz_1004.line ? xyz_1004 : xyz_1005
ind = get_ind(xyz,line,df_nav) # get Boolean indices
plot_path!(p1,xyz.traj,ind;show_plot=false,path_color=:black)
end
p1
end
Load previously processed maps
Eastern Ontario & NAMAD maps have been processed & saved previously. The full list of maps is in df_map.
begin # e_mask contains Boolean indices with "real" map data (not filled-in)
e_mapS_395 = get_map(:Eastern_395,df_map) # load map data
e_mask = trues(size(e_mapS_395.map)) # todo: update with v1.2.1
end;
begin
xx_lim = extrema(e_mapS_395.xx) .+ (-0.01,0.01)
yy_lim = extrema(e_mapS_395.yy) .+ (-0.01,0.01)
n_mapS_395 = upward_fft(map_trim(get_map(MagNav.namad),
xx_lim=xx_lim,yy_lim=yy_lim),e_mapS_395.alt)
end;
Plot all Ottawa area maps
3 maps are overlayed with 3 different color schemes.
begin
clims = (-500,500)
dpi = 50
p2 = plot_map(n_mapS_395;clims=clims,dpi=dpi,legend=false) # 395 m
plot_map!(p2,e_mapS_395;clims=clims,dpi=dpi,map_color=:magma) # 395 m
plot_map!(p2,p_mapS_800;clims=clims,dpi=dpi,map_color=:gray) # 800 m
end
Plot Eastern Ontario altitude map CDF
Most of the map data was collected at altitudes between 215 & 395 m.
begin # the Eastern drape map contains an additional drape (altitude) map
e_mapS_drp = get_map(:Eastern_drape,df_map) # load map data
e_alts = e_mapS_drp.alt[e_mask]
alt_avg = round(Int,mean(e_alts))
alt_med = round(Int,median(e_alts))
alt_val = 200:500
alt_cdf = [sum(e_alts .< a) for a in alt_val] / sum(e_mask)
p3 = plot(alt_val,alt_cdf,xlab="altitude [m]",ylab="fraction [-]",
title="altitude map CDF",lab=false,dpi=200)
end
Plot Eastern Ontario altitude map
Minimal areas have map data collected at 395 m or higher (colored spots), so a level map can be generated at 395 m using almost entirely upward continuation & minimal downward continuation.
begin
p4 = plot_map(e_mapS_395;dpi=dpi,legend=false,map_color=:gray)
e_mapS_395_ = deepcopy(e_mapS_395)
e_mapS_395_.map[e_mapS_drp.alt .<= 395] .= 0
plot_map!(p4,e_mapS_395_;dpi=dpi)
end
The trend of the map values agree, as expected.
Create a 3D map
All of the functionality in MagNav.jl can use 3D maps, which contain multiple stacked 2D maps with constant altitude spacing between each map level. There are 2 ways to create a 3D map, the first of which is to upward/downward continue to multiple altitudes using the upward_fft function. In this case, a filled Perth map at 800 m is stacked with the same map upward continued to 810 m.
begin
p_mapS3D_1 = upward_fft(p_mapS_800, p_mapS_800.alt .+ [0,10]) # 3D map
println(p_mapS3D_1.alt) # map altitude levels
p9 = plot_map(p_mapS3D_1) # plot single map level (default is lowest)
end
The second way to create a 3D map is to combine multiple 2D maps using the map_combine function. In this case, the Eastern map (395 m) & Perth map (800 m) are stacked into a 3D map.
begin
p_mapS3D_2 = map_combine([e_mapS_395,p_mapS_800]; # 3D map
n_levels = 2,
dx = MagNav.get_step(p_mapS_800.xx),
dy = MagNav.get_step(p_mapS_800.yy),
xx_lim = extrema(p_mapS_800.xx),
yy_lim = extrema(p_mapS_800.yy));
println(p_mapS3D_2.alt) # map altitude levels
p10 = plot_map(p_mapS3D_2) # plot single map level (default is lowest)
end